World Generation Experiments

These are my original (tm) experiments with terrain generation and erosion simulation. Of course I didn't invent those topics, but this is my take on them. I'll try and cite everywhere I get my crud from, but most of the topics are self citing. Like Perlin noise, the name is the citation.

Perlin Noise Baby!

from perlin_noise import PerlinNoise
p = PerlinNoise()
e = PerlinNoise(octaves=15)

p([0.9, 1])
e([0.9, 1])

-0.06844757119580991

Basic Rendering Setup

from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
import numpy as np
import random
%matplotlib widget
#plt.clf()
#plt.close()

def show_heightmap(z):
    x2, y2 = np.meshgrid(range(z.shape[0]), range(z.shape[1]))
    # show hight map in 3d
    fig = plt.figure()
    fig.set_figwidth(5)
    fig.set_figheight(5)
    ax = fig.add_subplot(projection='3d')
    ax.plot_surface(x2, y2, z)
    plt.show()

def show_image(z, cmap='grey'):
    fig = plt.figure()
    ax = fig.add_subplot()
    ax.imshow(z, cmap=cmap)
    plt.show()
    

def show_quiver(x_components, y_components):
    fig = plt.figure()
    ax = fig.add_subplot()
    ax.quiver(x_components, y_components)
    plt.show()

res = 300
det = 4

noise_generator = []
for i in range(det):
    noise_generator.append(PerlinNoise(octaves=((i+1)**2)))

zs = np.zeros([det, res, res])
for i in range(res):
    for j in range(res):
        for k in range(det):
            zs[k][i][j] += noise_generator[k]([i/res, j/res])
            #z[i][j] += noise_generator[k]([i/(k*res), j/(k*res)])/(k)
        

adding them normally makes it look fucked up

z = sum(zs)
show_heightmap(z)

Figure

so instead we add them by scaling each layer to be smaller and smaller. So the first octave represents big details, and thus has a big amplitude. By the time we get to the last octaves, we're making small details, so we give those layers smaller amplitudes.

This is a relatively simple idea, but its really important to terrain generation. It's actually this "fractal" way of adding things that makes stuff look good and natural, perlin noise itself is pretty bland, and there's a number of other noise algorithms that can make similar results- so long as we add them fractally.

z = np.zeros([res, res])
for i, hm in enumerate(zs):
    z += hm * (1.0/( (i+1)**2 ))
show_heightmap(z)

Figure

Realtime Procedural Terrain Generation (2004)

Defining Erosion

• Uses celluar automata to generate terrain, awesome.
• "Erosivness" is a score per each cell of the max difference between its height and its Von Neumann neighbors

grid_size = 256
cells = np.zeros([grid_size, grid_size])
def cell_erodedness(x, y):
    neighbors = np.array([
        cells[x][y+1],
        cells[x][y-1],
        cells[x+1][y],
        cells[x-1][y]
    ])
    return np.max(np.abs(neighbors - cells[x][y]))

• We want a map with lots of flat areas so that we can place game structures on them, but we don't want the map to be totally flat and thus boring.
• We'll achieve that by using the mean and mean $\bar{s}$ and variance $\sigma_s$
• Then we say the errosion score $\epsilon$ is $\epsilon = \frac{\sigma_s}{\bar{s}}$

def s_bar(cells_2d):
    return sum(cells_2d.flatten())/(grid_size**2)

def sigma_s(cells_2d):
    sb = s_bar(cells_2d)
    v = 0
    for cell_height in cells_2d.flatten():
        v += (cell_height - sb)**2
    return v / (grid_size**2)

def epsilon(cells_2d):
    return s_bar(cells_2d)/sigma_s(cells_2d)
    

Generating Base Terrain

• We want the base terrain to be close to erroded terrain to avoid having to run our errosion algorithm for too long (its computationally heavy).
• We'll start with pink noise (sometimes called $1/f$ noise)
• Then we'll add Voronoi diagrams and pertubation filtering to make it look more eroded.

Pink Noise

Pink noise also called $1/f$ noise is a fractal noise (just like Perlin) where the power spectral density is inversely porportional to the frequency.
And can be defined in 2 dimensions as this:

Or in engineer language, the amplitude of the higher-frequency sine waves is smaller than the lower frequency sine waves.

So a sine wave of frequency 2hz has an amplitude half that of a sine wave with frequency 1hz.

Pink noise appears a lot in nature, it is the sound waterfalls make, it appears as electrical signals in the human brain, it is similar to brownian motion which can be used to model the path of particle moving in a fluid. And of course, it can model the height map of mountains which have yet to be erroded.

Generating Pink Noise Via Fourier Transform

def gen_pink(size=128, scale=0.25, seed=42):
  # Generate white noise
  #rng = np.random.default_rng(int(seed))
  white = np.random.normal(size=(size, size))

  # Fourrier transform shift to centre lower frequency
  white_ft = np.fft.fftshift(np.fft.fft2(white))

  # Create a field such that as we go further from the origin
  # frequency goes up
  hsize = size / 2 # since we're centering [0, size] -> [-size/2, +size/2]
    
  # np.hypot finds the euclidean distance to the center (hypoteneuse)
  # np.ogrid just makes two vectors [-size/2, size], one vertical one
  # horizontal, the final result in a 2d matrix of distance to the middle
  # x[i][j] = distance (i,j) -> (0, 0)
  # uncomment the imshow below if you're still confused
  freq = np.hypot(*np.ogrid[-hsize:hsize, -hsize:hsize])
  #plt.imshow(freq)
  #plt.show()

  # Now we scale the noise in fourier space so that its frequency matches
  # the definition of pink noise (its squared to brownian noise maybe?)
  pink_ft = np.divide(white_ft, freq**2,
    out=np.zeros_like(white_ft), where=freq!=0)

  # Leave fourrier space
  pink = np.fft.ifft2(np.fft.ifftshift(pink_ft)).real

  # Scale the noise to be between 0 and 1
  lo, hi = np.min(pink), np.max(pink)
  pink = (pink - lo) / (hi - lo)
  return pink

show_image(gen_pink(), cmap='grey')

Figure

Using The Diamond-Square Algorithm

Pink noise can also be computed using the diamond square algorithm. The diamond square algorithm.

import math
import numpy as np
import random
def init_diamond_square(size):
    # Start with a terrain of height -1 everywhere
    terrain = np.zeros([size, size], dtype='float32') -1
    
    # Fill in the corners with uniform random values
    terrain[0][0] = np.random.random()
    terrain[0][size-1] = np.random.random()
    terrain[size-1][0] = np.random.random()
    terrain[size-1][size-1] = np.random.random()

    return terrain

def roughness_curve(roughness, avg, i, j):
    return (roughness * np.random.normal(1)) + (1.0 - roughness) * avg

def diamond_step(terrain, step_size, roughness):
    half_step = math.floor(step_size/2.0)
    x_steps = range(half_step, terrain.shape[0], step_size)
    y_steps = x_steps[:]
    for i in x_steps:
        for j in y_steps:
            if terrain[i][j] == -1:
                new_val = 0
                new_val += terrain[i-half_step][j-half_step]
                new_val += terrain[i-half_step][j+half_step]
                new_val += terrain[i+half_step][j-half_step]
                new_val += terrain[i+half_step][j+half_step]
                new_val /= 4
                terrain[i][j] = roughness_curve(roughness, new_val, i, j)
    return terrain

def square_avg(terrain, i, j, half_step):
    valid_neighbors = 4.0
    avg = 0.0

    if i - half_step >= 0:
        avg += terrain[i-half_step, j]
    else:
        valid_neighbors -= 1

    if i + half_step < terrain.shape[0]:
        avg += terrain[i+half_step, j]
    else:
        valid_neighbors -= 1

    if j - half_step >= 0:
        avg += terrain[i, j-half_step]
    else:
        valid_neighbors -= 1

    if j + half_step < terrain.shape[1]:
        avg += terrain[i, j+half_step]
    else:
        valid_neighbors -= 1

    return avg/valid_neighbors
        

def square_step(terrain, step_size, roughness):
    half_step = math.floor(step_size/2)
    steps_x_vert = range(half_step, terrain.shape[0], step_size)
    steps_y_vert = range(0, terrain.shape[1], step_size)
    steps_x_horiz = range(0, terrain.shape[0], step_size)
    steps_y_horiz = range(half_step, terrain.shape[1], step_size)

    #find_avg = lambda i, j: (terrain[i+1][j] + terrain[i-1][j] + terrain[i][j+1] + terrain[i][j-1])/4.0
    find_avg = lambda i, j: square_avg(terrain, i, j, half_step)
    
    for i in steps_x_horiz:
        for j in steps_y_horiz:
            terrain[i][j] =  roughness_curve(roughness, find_avg(i, j), i, j)
            
    for i in steps_x_vert:
        for j in steps_y_vert:
            terrain[i][j] = roughness_curve(roughness, find_avg(i, j), i, j)
            
    return terrain


def gen_diamond_square(x=8, roughness=1):
    size = 2**x + 1
    terrain = init_diamond_square(size)
    for i in range(x):
        # Smaller and smaller steps as we go
        step_size = math.floor( (terrain.shape[0]-1) / math.pow(2, i) )

        # Scale the roughness exponentially to match our increasingly smaller
        # diamonds and squares
        r = math.pow(roughness, i)
        
        # Once we run out of resolution quit
        if step_size == 0:
            break

        terrain = diamond_step(terrain, step_size, r)
        terrain = square_step(terrain, step_size, r)
    
    return terrain[:2**x,:2**x]
                              

a = gen_diamond_square(x=8, roughness=0.5)
show_image(a, cmap='grey')

Figure

# Lets run a test to see if any bad patterns emerge
s = np.zeros([2**8, 2**8])
sample_size = 50
for _ in range(sample_size):
    s += gen_diamond_square(x=8, roughness=0.5)
s /= sample_size
show_image(s)

Figure

show_heightmap(a)

Figure

• Although the Diamond-Square algorithm is much slower, it has the benefit of allowing us to change the random offset depending on the coordinates, so we can do things like simulate deposited material on less-steep terrain.
• Essentially we can re-write roughness_curve to include make use of the coordinates and create smoother or rougher areas. For example, below we great a mountain range that gets progressively smoother as approach the edges.

old_roughness_curve = roughness_curve
def roughness_curve(roughness, avg, i, j):
    roughness = roughness * 1.0/(np.hypot(i-64, j-64)+1)
    return (roughness * np.random.normal(1)) + (1.0 - roughness) * avg

a = gen_diamond_square(x=7, roughness=1)
show_heightmap(a)

Figure

• The above isn't really useful, but it demonstrates the effect we'll be using to create sediment from our supposed errosion.

# Set it back to something that isn't stupid lol
roughness_curve = old_roughness_curve

Voronoi Diagrams

• Vornoi Diagrams create cool cell-like natural structures, they're used in procedural textures all the time to create things like tissue, sponges, scales, pebbles, flagstones, and for our case, mountains.

• They're actually really simple, first we create $N$ features points, which are given a random x and y. Then we calculate the 2D euclidean distance from each cell to each point, and take some weighted sum of them.

$$ h = c_1d_1 + c_2d_2 + c_3d_3 + ... c_nd_n $$

• After fiddling I found out that its also common to just take the distance to the maximum value, which produces a very similar result and runs much faster.

def gen_voronoi(size, n_points):
    points_x = np.random.uniform(0, size, size=n_points)
    points_y = np.random.uniform(0, size, size=n_points)
    c1 = -1 # special to this paper
    c2 = 1
    cells = np.zeros([size, size])
    for i in range(size):
        for j in range(size):
            distances = ((points_x - i)**2 + (points_y - j)**2)
            # First closest gets multiplied by coefficient 1
            min1 = np.argmin(distances)
            cells[i][j] = min(distances)*c1
            # Then make him really big so we can find the second distance
            distances[min1] = size**2
            # and multiply it by coefficient 2
            cells[i][j] += min(distances)*c2

    # Rescale to be between 0 and 1
    cells -= min(cells.flatten())
    cells /= max(cells.flatten())
    return cells

a = gen_voronoi(256, 64)
show_image(a)

Figure

Combination and Perturbation

• We just came up with a bunch of really nice looking noise functions. Now we get to combine them to make cool stuff!
• We'll combine 2/3rds diamond-square with 1/3rd voronoi diagram.
• Then we'll use "bomb" style random placement patterns, which don't come from this paper, but a book it references.

def gen_combined(size, ratio=0.33, roughness=0.25, peaks=16):
    ds = gen_diamond_square(x=int(math.log(size, 2)), roughness=roughness)*(1-ratio)
    vo = gen_voronoi(size, peaks )*ratio
    return ds + vo

from ipywidgets import interact, interactive, fixed, interact_manual
def interact_gen_combined(ratio=0.5, roughness=0.5, height=1.0):
    a = gen_combined(512, ratio=ratio, roughness=roughness)
    a = a*height
    show_image(a)
    show_heightmap(a)
interact_manual(interact_gen_combined)

interactive(children=(FloatSlider(value=0.5, description='ratio', max=1.5, min=-0.5), FloatSlider(value=0.5, d…

<function __main__.interact_gen_combined(ratio=0.5, roughness=0.5, height=1.0)>

Gradients and Quivers

Because of the way we've organized our heightmap we are unable to use np.gradient to compute the three dimensional gradient, but we can still compute it in 2 dimensions. It's sort of like having a compass that points to the nearest hole.

gr = np.gradient(a)

a = np.meshgrid(np.arange(0, 20), np.arange(0, 20))
a = np.random.normal(0, 1, size=[20, 20])
b = np.gradient(a)

fig = plt.figure()
ax = fig.add_subplot()
ax.quiver(b[0], b[1])
plt.show()

Figure

v = gen_diamond_square(x=7, roughness=0.5)
v.shape

(128, 128)

vg = np.gradient(v)

fig = plt.figure()
ax = fig.add_subplot()
ax.quiver(vg[0], vg[1])
plt.show()

Figure

dd = np.hypot(vg[0], vg[1])

show_image(dd)

Figure

from PIL import Image
im = Image.fromarray((dd*4096).astype(np.uint8), 'L')
#im.save('test.jpg')

Performance Concerns Regarding np.gradient

By default numpy's gradient function computes the gradient for the entire array. This is problametic because we're going to be re-computing the gradient almost every step, and we only really care about the gradient at 1 or 2 points throughout. To remedy this we'll look under the hood at what numpy is doing, and try to replicate it on a cell-by-cell basis.

import numpy as np

test_grid = np.random.normal(size=[1024, 1024])

import time
start = time.time()
for _ in range(500):
    np.gradient(test_grid)
print(time.time() - start)

5.394184350967407

def single_gradient(arr, i ,j):
    x_comp = (arr[i+1][j] - arr[i-1][j])/2.0
    y_comp = (arr[i][j+1] - arr[i][j-1])/2.0
    return x_comp, y_comp

example = np.random.normal(size=[128, 128])
single_gradient(example, 1, 1)

(0.38082777585925565, 1.235753511308505)

example3 = np.random.normal(size=[12, 12])
gr3 = np.gradient(example3)

single_gradient(example3, 3, 6)

(0.11610482642876868, -0.4033757916001024)

gr3[0][3][6], gr3[1][3][6]

(0.11610482642876868, -0.4033757916001024)

start = time.time()
for _ in range(500):
    single_gradient(test_grid, 2, 2)
print(time.time() - start)

0.0003571510314941406

z1 = np.zeros([128, 128]) z2 = np.zeros([128, 128]) for i in range(128): for j in range(128): z1[i][j], z2[i][j] = single_gradient(example, i, j)

Alright, that seems to have done the trick. Trust me, this will save us a lot of computing time later when we implement that actual erosion algorithm.

Particle Based Erosion

We're going to model drops of water falling onto the heightmap and running downhill. Instead of implementing 3D water droplets and having them move up and down a 3D heightmap, we'll have the particles operate in 2 dimensions, and we'll treat the heightmap like an acceleration field. So particles will accelerate in $X$ and $Y$ as it passes over a "slope", which will be given by the gradient of our heightmap. This is very similar to the way we sometimes model charged particles passing through magnetic fields.

Finding Downhill

Of course, to accelerate it downhill we'll need to know where "downhill" is. We'll do that using the gradient of the previous heightmap. To start, instead of using some noise function, let's make a simple radial hill, so that mistakes will be obvious.

def radial_hills(size):
    g = np.zeros([size, size])
    for i in range(size):
        for j in range(size):
            g[i][j] = ((i - size/2.0) ** 2 + (j - size/4.0) ** 2)**0.5
    return g

grid = radial_hills(20)

show_image(grid)
grad = np.gradient(grid)
show_quiver(-grad[1], -grad[0])

Figure

Figure

Actually Doing Erosion

For particle based erosion we will simualte drops of water falling on the terrian and moving the dirt (sediment) around. To do this we will use the chemical engineering concept of mass transfer coefficients. This will model the amount of dirt that is "in solution" vs "deposited".

Mass transfer is proportional to the difference between the concentration $c$ and the eqiilibirum concentration $c_{eq}$. $$ \frac{dc}{dt} = k(c_{eq} - C) $$ Where $k$ is the mass transfer coefficient. This equation means that the further we get from the equilibrium point, the faster we will deposite/dissolve sediment to reach it.The difference between the equilibrium point and the actual concentration is often reffered to as the driving force.

These particles will follow two dimensional newtonian physics. We'll spawn it randomly on the heightmap, as it bumps into steeper slopes, it will accelerate down the slope.

class Drop:
    def __init__(self, x, y, evaporation_rate=0.01, friction=0.01,  deposition_rate=0.5):
        self.pos = np.array([x, y], dtype='float64')
        self.speed  = np.array([0.0, 0.0])
        self.volume = 1.0 # how much water
        self.sediment = 0.0 # how much dirt
        self.inv_friction = 1.0 - friction # how much velocity to keep each tick
        self.evaporation_rate = evaporation_rate # how much water to keep each tick
        self.deposition_rate = deposition_rate # linearly scales dc/dt
        self.prev_height = None # set once the particle does its first erosion
        
    # Moves the particel following 2D newtonian motion
    def move(self, dt, accel):
        self.speed += dt*accel/(self.volume*0.1)
        self.pos += dt*self.speed
        self.speed *= self.inv_friction

    # Returns how much sediment to take off the heightmap
    def erode(self, dt, height_at_pos):
        # This stops the water from having an unaturally high effect on the
        # first spot it lands (because it starts with 0 dirt)
        if self.prev_height is None:
            self.prev_height = height_at_pos
            return 0.0

        # Compute the mass transfer coefficient at equilibrium
        c_eq = self.volume*np.hypot(*self.speed)*(height_at_pos - self.prev_height)
        if c_eq < 0.0:
            c_eq = 0.0

        # Then compute the "driving force", the difference between equilibrium and where we are
        cdiff = c_eq - self.sediment

        # Deposite a bit of sediment 
        self.sediment += dt*cdiff*self.deposition_rate
        
        # Evaporate a bit of water
        self.volume *= self.volume * (1.0 - self.evaporation_rate*dt)

        # Save this so we don't have to look it up again
        self.prev_height = height_at_pos

        # Return the difference in height that should have been made at that
        # spot on the terrain
        return dt*self.volume*self.deposition_rate*cdiff

    def should_kill(self, width, length):
        return self.pos[0] > width or self.pos[1] > length or self.pos[0] < 0 or self.pos[1] < 0 or self.volume < 0.01

def erode_heightmap(heightmap, dt=1.0, cycles=100, evaporation_rate=0.01, friction=0.01, deposition_rate=0.5):
    # To save some typing
    width, length = heightmap.shape
    dt = dt # how much time passes every tick
    
    # Make a random drop to begin with
    drop = Drop( 
            random.randint(0, width-1), random.randint(0, length-1),
            evaporation_rate = evaporation_rate,
            friction = friction,
            deposition_rate = deposition_rate
    )

    # Then for every "cycle"
    for cycle in range(cycles):
        # Round down the position so we can use it as the index of the gradient
        posx, posy = drop.pos
        posx, posy = int(posx -1), int(posy -1)

        # Add or remove sediment according to the drops equilibrium point
        heightmap[posx, posy] -= drop.erode(dt, heightmap[posx, posy])

        # Accelerate the drop along the gradient of it's current position
        accl = np.array(single_gradient(heightmap, posx, posy))
        
        if np.isnan(accl[0]) or np.isnan(accl[1]):
            raise ValueError("NaN in acceleration. Somebody's violating the conservation of matter and energy." + str( (posx, posy)))
            
        drop.move(dt, accl)

         # If the drop runs out of water or goes out of bounds, replace it with a new drop
        if drop.should_kill(width, length):
            drop = Drop(
                random.randint(0, width-1), random.randint(0, length-1),
                evaporation_rate = evaporation_rate,
                friction = friction,
                deposition_rate = deposition_rate
            )

gr = gen_pink(size=128)
og = np.copy(gr)
i=3
cycles = 1000000*6
erode_heightmap(gr, dt=0.25, cycles=cycles, evaporation_rate=0.00001, friction=0.0001, deposition_rate= (1.0/6.0)*(i+1))

f, axarr = plt.subplots(2, 2, figsize=[10, 10])
axarr[0][0].imshow(og, cmap='gray')
axarr[1][0].imshow(og - gr, cmap='gray')
axarr[0][1].imshow(gr, cmap='gray')
plt.show()

Figure

from ipywidgets import IntText, FloatSlider, IntSlider, FloatText
import time
grid = gen_pink(size=64)
original_grid = np.copy(grid)
cycle_count = 0

def interact_erode_heightmap(size=256, dt=1.0, cycles=1000, evaporation_rate=0.01, friction=0.01, deposition_rate=0.5):
    global grid, original_grid, cycle_count

    # If we change the size we'll have to reset the graph
    if grid.shape[0] != size:
        grid = gen_pink(size=size)
        original_grid = np.copy(grid)

    cycle_count += cycles
    print(size, dt, cycles, evaporation_rate, friction, deposition_rate, cycle_count)

    # Copy the original grid and run erosion on it
    start_gen_time = time.time()
    erode_heightmap(grid, dt=dt, cycles=cycles, evaporation_rate=evaporation_rate, friction=friction, deposition_rate=deposition_rate)
    # Print how long that took
    print("Generation time in seconds:", time.time() - start_gen_time)
    
    # Make a plot showing the original graph, eroded graph, the difference between them, and the gradient of the original
    f, axarr = plt.subplots(2, 2, figsize=[7, 7])
    axarr[0][0].imshow(original_grid, cmap='gray')
    axarr[1][0].imshow(original_grid - grid, cmap='gray')
    axarr[0][1].imshow(grid, cmap='gray')
    g = np.gradient(original_grid)
    axarr[1][1].quiver(g[1], g[0])
    plt.show()

interact_manual(
    interact_erode_heightmap,
    size=IntText(value=128),
    dt=FloatText(value=0.25),
    cycles=IntText(value=1000000*5),
    evaporation_rate=FloatText(value=0.00001),
    friction=FloatText(value=0.0001),
    deposition_rate=FloatText(value=0.66)
)

interactive(children=(IntText(value=128, description='size'), FloatText(value=0.25, description='dt'), IntText…

<function __main__.interact_erode_heightmap(size=256, dt=1.0, cycles=1000, evaporation_rate=0.01, friction=0.01, deposition_rate=0.5)>

Parameter Tuning

I'm not really sure what affect the deposition_rate has on the simulation, or really what the overall best paramater choices are. So let's try out a few and see what they look like.

f, axarr = plt.subplots(3, 2, figsize=[8, 8])
cycles = 100000*3
og_go = gen_pink(size=64)
go = np.copy(og_go)
for i in range(2*3):
    gr = np.copy(go)
    erode_heightmap(gr, dt=0.25, cycles=cycles, evaporation_rate=0.00001, friction=0.0001, deposition_rate= (1.0/6.0)*(i+1))
    axarr[i%3][i%2].imshow(gr, cmap='gray')
    axarr[i%3][i%2].text(0, 5, i, color='yellow')
plt.show()

Figure

Higher deposition rates seem to lead to more of those "streaky" channels, and also some tessolation- although that could just be from the low resolution. I like the way 1 looks, but they all seem acceptable.

Now let's try friction!

%store
f, axarr = plt.subplots(3, 2, figsize=[8, 8])
go = np.copy(og_go)
for i in range(2*3):
    gr = np.copy(go)
    erode_heightmap(gr, dt=0.25, cycles=cycles, evaporation_rate=0.0001, friction=1*(10**(-2-i)), deposition_rate=0.5)
    axarr[i%3][i%2].imshow(gr, cmap='gray')
    axarr[i%3][i%2].text(0, 5, i, color='yellow')
plt.show()

Stored variables and their in-db values:

Figure

The difference here is pretty hard to see to be honest.

Now let's try fiddling with the evaporation rate. I have a feeling it will do something similar to friction.

f, axarr = plt.subplots(3, 2, figsize=[8, 8])
go = np.copy(og_go)
for i in range(2*3):
    gr = np.copy(go)
    erode_heightmap(gr, dt=0.25, cycles=cycles, friction=(10**(-5)), evaporation_rate=1*(10**(-2-i)), deposition_rate=0.5)
    axarr[i%3][i%2].imshow(gr, cmap='gray')
    axarr[i%3][i%2].text(0, 5, i, color='yellow')
plt.show()

Figure

It looks like it does something similar to the deposition rate, except it doesn't seem to have the same tessolation problem. I like the way 5 looks (the lowest evaporation rate, and thus longest lasting drops). It seems to create more "erosiony" channels.

Optimizing by Combining with Previous Techniques

Those results look pretty nice (I rendered them in blender and they actually look very nice) but they have the same flaw as most erosion simulators, which is that erosion simulation is very expensive. Generating a terrain of 256x256 units alone takes a while- and doesn't give that much detail. Naturally we'll want larger more complete terrains, to do that we'll try and start with more eroded terrain than normal.

See how this already has similar patterns to our eroded terrain? Let's see if we can avoid running as many cycles on it.

show_image(gen_combined(128, ratio=0.4, roughness=0.47, peaks=16))

Figure

# Generate some with pure pink noise and some with our combined noise algorithm grp = gen_pink(size=128) gr = gen_combined(128, ratio=0.4, roughness=0.47, peaks=16*2) og = np.copy(gr) ogp = np.copy(grp) # Perform erosion on both i=3 cycles = (1000000*6)//2 erode_heightmap(gr, dt=0.25, cycles=cycles, evaporation_rate=0.00001, friction=0.0001, deposition_rate= (1.0/6.0)*(i+1)) erode_heightmap(grp, dt=0.25, cycles=cycles, evaporation_rate=0.00001, friction=0.0001, deposition_rate= (1.0/6.0)*(i+1)) # Draw them so we can compare f, axarr = plt.subplots(2, 2, figsize=[10, 10]) # Top combined axarr[0][0].imshow(og, cmap='gray') axarr[0][1].imshow(gr, cmap='gray') # Bottom pnk axarr[1][0].imshow(ogp, cmap='gray') axarr[1][1].imshow(grp, cmap='gray') plt.show()im = Image.fromarray((grp*200).astype(np.uint8), 'L') im.save('heightmap.png')

Yeah the pure pink looks WAAAY nicer. Maybe that's just because my Diamond-Square pink noise generator is kind of shit?

# Generate some with pure pink noise and some with our combined noise algorithm grp = gen_pink(size=128) gr = 0.5*gen_voronoi(128, 16) + 0.5*grp og = np.copy(gr) ogp = np.copy(grp) # Perform erosion on both i=3 cycles = (1000000*6)//2 erode_heightmap(gr, dt=0.25, cycles=cycles, evaporation_rate=0.00001, friction=0.0001, deposition_rate= (1.0/6.0)*(i+1)) erode_heightmap(grp, dt=0.25, cycles=cycles, evaporation_rate=0.00001, friction=0.0001, deposition_rate= (1.0/6.0)*(i+1)) # Draw them so we can compare f, axarr = plt.subplots(2, 2, figsize=[10, 10]) # Top combined axarr[0][0].imshow(og, cmap='gray') axarr[0][1].imshow(gr, cmap='gray') # Bottom pnk axarr[1][0].imshow(ogp, cmap='gray') axarr[1][1].imshow(grp, cmap='gray') plt.show()

Pseudo-Chunking

By Linear Interpolation:

def make_eroded_chunk(size):
    grp = gen_pink(size=size)
    i=3
    cycles = (1000000*6)//4
    erode_heightmap(grp, dt=0.25, cycles=cycles, evaporation_rate=0.00001, friction=0.0001, deposition_rate= (1.0/6.0)*(i+1))
    return grp

n = 64 # size of the chunk
m = 8  # size of the seam

chunk_a = make_eroded_chunk(n + 2*m)
chunk_b = make_eroded_chunk(n + 2*m)

1. Generate a chunk of size $n+2m$, call it chunk (0, 0). Erode it for $E$ cycles.
2. Generate a second chunk of size $n+2m$, call it chunk (1, 0). Erode it for $E/2$ cycles.
3. Pad the right of chunk (0, 0) with 0s, until it $2*n + 2*m$ in the X direction.
4. Pad the left of chunk (1, 0) with 0s similarly.
5. Add chunk (0, 0) and chunk (1, 0) together. Call this the chunk quilt.
6. At the overlap site $x=[n, n+2m]$, linear interpolate between the two chunks.
7. Erode the from position $x=[n, 2m+2m]$ only.

First, we'll try the obvious way of blending these two shapes together by linear intropolating them.

left_chunk = np.concatenate( (chunk_a, np.zeros([n+2*m, n])), axis=1)
right_chunk = np.concatenate( (np.zeros([n+2*m, n]), chunk_b), axis=1)
chunk_quilt = left_chunk + right_chunk

mask_a = np.ones(chunk_quilt.shape)
mask_b = np.ones(chunk_quilt.shape)
for i in range(n, n+2*m):
    for j in range(0, n+2*m):
        mask_a[j, i] = (i-n)/(2*m)
        mask_b[j, i] = 1.0 - ((i-n)/(2*m))
show_image(left_chunk*mask_b + right_chunk*mask_a)

Figure

That honestly doesn't look half bad, but we can clearly see some blurry staircase-like patterns at the seam, that's no good! The fundamental problem with this approach is that the two chunks dont have anything naturally in common, so it'll take a lot of work to force one into the other, and it'll always look somewhat unnatural.

By progressive erosion

This technique is based off the theory that pink noise cells generated next to each other will bare some resemblance (or cooperation?) two eachother, even if eroded separately. We'll call these adjacent pink noise chunks "sister chunks" because it sounds cool.

# Let's agree on some erosion parameters
drops_per_pixel_constant = 2000000 / ((2*64 + 2*8) * (64 + 2*8)) # from a sample n=64 n=8 that looked good
cycles = int(drops_per_pixel_constant*((2*n+2*m) * (n+2*m)))

# Erosion settings
quick_erode = lambda c: erode_heightmap(
    c, dt=0.25, cycles=cycles, evaporation_rate=0.00001, friction=0.0001, deposition_rate=0.66
)

# First let's try it without any "seam buffer" (m=0)
mega_chunk = gen_pink(size=(2*n ))
mega_chunk = mega_chunk[0:n, 0:]
left_chunk = mega_chunk[0:, 0:n]
right_chunk = mega_chunk[0:, n:]
show_image(np.concatenate( (left_chunk, right_chunk), axis=1))
quick_erode(left_chunk)
quick_erode(right_chunk)
show_image(np.concatenate( (left_chunk, right_chunk), axis=1))

Figure

/tmp/ipykernel_38863/1435930297.py:2: RuntimeWarning: More than 20 figures have been opened. Figures created through the pyplot interface (`matplotlib.pyplot.figure`) are retained until explicitly closed and may consume too much memory. (To control this warning, see the rcParam `figure.max_open_warning`). Consider using `matplotlib.pyplot.close()`.
  fig = plt.figure()

Figure

Obviously this isn't an acceptable end-product, but it does confirm the theory about sister chunks eroding similarly, even when done separately. They're way more similar than I expected, and it might be possible simply to stitch these together with linear interpolation.

By linear interpolating progressive erosion

# Chunk size and seam size
n = 64
m = 32

# We'll try to guess how much rain an area of this size should get
drops_per_pixel_constant = 3000000 / ((2*64 + 2*8) * (64 + 2*8)) # from a sample n=64 n=8 that looked good
cycles = int(drops_per_pixel_constant*((2*n+2*m) * (n+2*m)))

# Erosion settings
quick_erode = lambda c: erode_heightmap(
    c, dt=0.25, cycles=cycles, evaporation_rate=0.00001, friction=0.0001, deposition_rate=0.66
)

# Generate two pink noise chunks together
mega_chunk = gen_pink(size=(2*n + 2*m))
mega_chunk = mega_chunk[0:n+2*m, 0:]

# Left chunk and right chunk will share a seam buffer
left_chunk = np.copy(mega_chunk[0:, 0:n+2*m])
right_chunk = np.copy(mega_chunk[0:, n:])

# Erode them individually
quick_erode(left_chunk)
quick_erode(right_chunk)

# Pad each one to the size of the whole terrain (this is inefficient but easy to work with)
pad_left = np.concatenate( (left_chunk, np.zeros([n+2*m, n])), axis=1)
pad_right = np.concatenate((np.zeros([n+2*m, n]), right_chunk), axis=1)

# Linear interpolate between them
mask_a = np.ones([n+2*m, 2*n + 2*m])
mask_b = np.ones([n+2*m, 2*n + 2*m])
for i in range(n, n+2*m):
    for j in range(0, n+2*m):
        mask_a[j, i] = (i-n)/(2*m)
        mask_b[j, i] = 1.0 - ((i-n)/(2*m))

# Show off
show_image(pad_left*mask_b + pad_right*mask_a)

Figure

show_image(pad_left*mask_b + pad_right*mask_a)

Figure

# Show where the seams are
plt.axvline(x=n+m, color='red', ls='--')
plt.axvline(x=n, color='yellow', ls='--')
plt.axvline(x=n+2*m, color='yellow', ls='--')

<matplotlib.lines.Line2D at 0x7751a2b80ec0>

Predetermined ridges and Voronoi Problem

This idea here was to pick the Voronoi points for the whole world (or pseudorandomly generate them, then generate the ones close enough to matter), then presume they will become the ridges of the mountain. So when chunk A generates, it knows where the ridges of chunk B are, even if chunk B hasn't even been generated yet.

I've decided not to pursue this one, because it's a more-deliberate version of the pink noise solution above. In the above, the ridges are determined by the location of the high-points in the pink noise, and while they don't always generate in exactly the same spot, it seems to work well enough to hide the seam.

If however, we decide to add back the Voronoi noise, it may create a problem with those ridges being unkown (because those Voronoi points haven't been chosen yet).

So there's a potential problem with Voronoi noise being un-chunked, but seeding the random number generation of the points should allow us to re-chunk it again.

Final Implementation

Of course, Pink Noise being based off a lot of trig functions, is recurring, and not very good for an "endless chunked map". We can get around this by just generating a really big pink noise map and saving it, then "revealing" it as we go- but that sort of defeats the "infinite terrain" fun of a chunking algorithm.

Instead, I'll just swap it out for Perlin Noise, they look somewhat similar right?

from perlin_noise import PerlinNoise
def perlin_chunk(x=0, y=0, size=64, seed=42, octaves=8):
    x, y = float(y)*size, float(x)*size
    waves = [ PerlinNoise(octaves=i+1, seed=seed+i) for i in range(octaves) ]
    #waves = [ PerlinNoise(octaves=octaves, seed=seed) ]
    chunk = np.zeros( [size, size] )
    for i in range(size):
        for j in range(size):
            chunk[i, j] = sum([ w( (((i+x)/size), ((j+y)/size)) ) for w in waves ])
    return chunk

show_image(np.concatenate([perlin_chunk(), perlin_chunk(x=1)], axis=1))

Figure

Then we'll write a function that generate two "sister chunks" from the Perlin noise, and erodes them in similar ways:

def gen_eroded_chunk(x=0, y=0, size=64, seed=42, octaves=8):
    # Generate the perlin noise chunk at this position
    c = perlin_chunk(x=x, y=y, seed=seed, octaves=8, size=size)
    
    # We'll try to guess how much rain an area of this size should get
    drops_per_pixel_constant = 3000000 / ((2*64 + 2*8) * (64 + 2*8)) # from a sample n=64 n=8 that looked good
    cycles = int(drops_per_pixel_constant*(size**2))
    
    # Erosion settings
    quick_erode = lambda c: erode_heightmap(
        c, dt=0.25, cycles=cycles, evaporation_rate=0.00001, friction=0.0001, deposition_rate=0.66
    )

    # Then just perform simple erosion on it
    quick_erode(c)
    return c

a, b = gen_eroded_chunk(), gen_eroded_chunk(x=1)
show_image(np.concatenate( [a, b], axis=1))

Figure

Just like before, the major features are okie-dokie, but the minor ones make a clear seam. We'll hide this with a linearly-interpolated seam buffer like before. I'll write a function to stitch the right and left chunks, and then we can just rotate/transpose the matrices to make them fit enough to stitch the other chunks.

def stitch_chunk_right(chunka, chunkb, seam_size=8):
    seam = np.zeros( [chunka.shape[0],  seam_size] )
    seamx_start = chunka.shape[0] - seam_size
    for i in range(seam_size):
        for j in range(chunka.shape[1]):
            seam[j, i] = chunka[j, seamx_start+i]*(1.0 - i/seam_size) + chunkb[j, i] * (i/seam_size)
    return chunka[0:chunka.shape[1], 0:seamx_start], seam, chunkb[0:chunkb.shape[1], seam_size:]

n_a, seam, n_b = stitch_chunk_right(a, b, seam_size=32)
quilt = np.concatenate( (n_a, seam, n_b), axis=1)
show_image(quilt)
print(quilt.shape)

Figure

(64, 96)